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Abstract. This article presents a new mathematical framework to perform statis- 
tical analysis on time-indexed sequences of 2D or 3D shapes. At the core of this 
statistical analysis is the task of time interpolation of such data. Current models in 
use can be compared to linear interpolation for one dimensional data. We develop a 
spline interpolation method which is directly related to cubic splines on a Riemannian 
manifold. Our strategy consists of introducing a control variable on the Hamiltonian 
equations of the geodesies. Motivated by statistical modeling of spatiotemporal data, 
we also design a stochastic model to deal with random shape evolutions. This model 
is closely related to the spline model since the control variable previously introduced 
is set as a random force perturbing the evolution. 

Although we focus on the finite dimensional case of landmarks, our models can 
be extended to infinite dimensional shape spaces, and they provide a first step for a 
non parametric growth model for shapes taking advantage of the widely developed 
framework of large deformations by diffeomorphisms. 

1. Introduction 

Mathematical and statistical modeling of shapes has undergone significant development 
over the past twenty years driven by a wide range of applications in medical imaging. Ini- 
tially the focus was on the comparison between two shapes also refereed to as registration. 
Among others, registration and comparison tools derived from a Riemannian point of view 
on shapes spaces and diffeomorphic transport have been actively developed during the past 
few years. This framework was used to represent shapes and study the statistical vari- 
ation of static shapes within a population. An emerging question of interest is now to 
study time dependent data of shapes (images, landmarks, surfaces or tensors). The basic 
dataset is then a sequence of shapes indexed by time. For example, a practical application 
would be the analysis of follow-up studies in brain imaging. 

Several attempts models for the variability of longitudinal data have been proposed 
recently: a parametric model of growth is proposed in [Ml [15], which aims to describe 
the biological evolution as an iteration of random elementary diffeomorphisms, so called 
GRID. Focusing on image data, statistical estimation of the parameters are performed 
with the GRID model in [3Q| 132] , Other attempts are often based on using an initial 
registration tool (such as geodesies on a group of diffeomorphisms, see [35] for a large 
overview) to interpolate the time dependent data with piecewise geodesies [9l|22]. In fO^, 
the model is further developed with the introduction of time realignments to allow the 
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study of an ensemble of longitudinal data and the computation of an averaged space-time 
evolution. 

From the modeling point of view a typical growth evolution is usually smooth in time. 
However, the piecewise geodesic models underlying current analysis of time dependent 
shape data can not prevent a loss of regularity at the observation points. Moreover, from 
a more classical statistical point of view, piecewise linear regression does not provide the 
best interpolation framework and from a probabilistic point of view the limiting process 
underlying piecewise geodesic interpolation is more or less a Kunita flow [53] that has a 
Brownian like time evolution. These remarks suggest that a better model of interpolation 
in time should be of higher order than piecewise geodesies, and also closely related to a 
probabilistic model of random evolution of shapes. The question of time interpolation 
was addressed in [TJ with the use of a kernel on the time variable but still with underlying 
piecewise geodesies on the data. Random evolutions of shapes have been treated for 
instance in [HI [21], but here again the model is of first order in the sense that the evolution 
is not smooth in time. 

In this work, we present a second order model for deterministic and stochastic shape 
evolution as a primary step toward further statistical applications. This work will ex- 
tensively use the Hamiltonian framework that emerged several years ago in the field of 
Computational Anatomy (CA) to compute the registration of landmarks (points of in- 
terest on an image) [T3l \T9\ 126] , This problem of landmark registration is equivalent to 
finding geodesies on the Riemannian manifold of landmarks. Thus shooting methods as 
well as gradient descent have been applied to solve this problem [4^,128]. 

A first natural step to deal with random shape evolution is to study a stochastic per- 
turbation of the Hamiltonian equations of geodesies. If we consider the evolution of 
landmarks as a physical system of particles, it corresponds to the introduction of a ran- 
dom force on their evolution providing smooth time random perturbation around a mean 
geodesic trajectory. However, this in turn leads to a new deterministic counterpart when 
the stochastic term in the evolution of the momentum is replaced by a deterministic con- 
trol variable. Given a sequence of time-indexed shapes, and optimizing on the control 
variable, we end up with a framework for smooth time interpolation between shapes. In 
the finite dimensional case of landmarks, this approach is directly related to splines on 
a Riemannian manifold [HI 123 UHl E] • We will present a Hamiltonian approach to derive 
the equations for the splines in a Hamiltonian setting. It can be related to the work of 
[TT] since they use a control point of view to obtain the Hamiltonian equations but our 
approach differs from theirs in that we use the Hamiltonian formulation of geodesies on 
landmark space. Hence there is a straightforward generalization to other Hamiltonian 
systems. The possible extension of the Hamiltonian equations of geodesies to curves, sur- 
faces (or embeddings of manifolds in R^) is the key feature to design second order models 
of evolution on continuous shape spaces [131 ESI [34] but in this work we will concentrate 
on the finite dimensional case of landmarks. Note however, that contrary to the initial 
context of air-craft trajectory planning (where the target manifold is the low-dimensional 
Lie Group SE(3) - the Special Euclidean group on R^), for which splines on Riemannian 
manifold were introduced - shape spaces in Computational Anatomy are more strongly 
connected with the infinite dimensional case. Moreover, we will consider both the deter- 
ministic and the stochastic setting since these are considered as fundamental ingredients 
towards statistical analysis of longitudinal data. 

The plan of the paper is as follows: we first introduce in section |2] the Hamiltonian 
equations for the case of landmarks, which constitute the central tool in this article. 
The section [3] is devoted to the spline model with different metrics and section [4] gathers 
theoretical results about the global existence in time of the second order deterministic 
spline evolution, the existence of a minimizer for the spline estimation problem and the 
derivation of the Euler-Lagrange equations. Then in section [5j we present some numerical 
experiments about shape spline estimation on 2D shape evolution, highlighting some of 
their main features. In section [6j we turn to problem of stochastic second order spline 
evolutions with a proof of the well-posedness and global existence in time and a few 
illustrative simulations. Last we develop in section [7| a slightly more general picture of 
shape splines on homogeneous spaces extending the model beyond the landmark setting. 
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2. Hamiltonian equations of geodesics for landmark matching 

The problem of landmark matching via diffeomorphic transport is now quite well un- 
derstood. The basic idea is to build a continuous path of minimal length (j)t starting 
from 00 = Idj^d in a group Gy of diffeomorphisms, between an initial configuration 
X = (^i)i<i<n of ^ landmarks in and a target configuration y = {yi)i<i<n • Thus, if 
Xi^t = and Xt = {xi^t)i<i<m is a path from x to y in the space of landmark con- 

figurations induced by a geodesic in a Riemannian space of diffeomorphisms. The group 
Gy is defined through the fiow of time dependent velocity fields (t, x) Vt{x) on 

(1) a^ = "*°* 

where F is a Hilbert space of velocity fields and v G L^([0,l],y) is an element of the 
space of time dependent velocity fields with finite norm. Obviously, the existence of a 
fiow t ^ for V e i^^([0, 1], V) solution of ^ depends on some regularity assumptions 
of the instantaneous velocity fields Vt, mainly the control of the first order derivatives 
of Vt. Assuming this, the group Gy is defined as Gy = {^i | v ^ L^([0, and 
the diffeomorphic matching problem for landmarks is formulated through the following 
variational problem 



min J^\vt\^ydt, veL\[0,l],V) 

(2) { with 

^\{xi) = yi, l<i<n. 

The problem ^ is well posed as soon as V is continuously embedded in Cq (M^, M^) the 
space of velocities vanishing at oo (such a V is called an admissible space). Namely, 
we assume the existence of a constant C such that for any v G V, the following inequality 
is verified 

(3) |^|l,oo <C|^|y. 

Such an admissible space is a Reproducible Kernel Hilbert Space (RHKS) and is equipped 
with a kernel Ky{z, z') playing a key role in defining the structure of the solution of the 
geodesic emerging from Since the reader may not be familiar with RKHS, let us say in 
a nutshell that for any G the kernel Ky{z^ z') is a d x d matrix in M(^(R), that z 
Ky{z,z')a' e V for any z',a' G R^, that {Ky{., z')a\ Ky{., z'')a'')y = a'^ Ky{z' , z")a" 
(the so called reproducible property) and that {v^Ky{.^ z')a')y = {v{z')^a')y for any 
V . The main fact is that given Ky^ the space V is completely defined and one may 
start from the choice of the kernel Ky itself to define the space V. A large number of 
kernels have been proposed most commonly the Gaussian kernel 

(4) Kv{z,z') = exp(- '^~f '' )IdR. . 

Even more importantly, the kernel appears explicitly in geodesics emerging from Q as 
described by the following Theorem. To ensure the existence of a solution in this theorem, 
we need to assume that the kernel is positive definite: 



(5) I ^ Ky{xi, .)ai\l = ^ V z = . 
Theorem 1 (cf [13j). The solution v G i^^([0, 1], V) exists and satisfie 

n 

(6) Vt{z) = Ky{z, Xi^t)Pi,t 

i=l 

where t ^ {xt^Pt) is solution of 
(7) 

with Ho{p, x) = j pfKvixi, Xj)pj and xo = x. 
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An interesting fact is that the evolution ([7|) is a Hamiltonian evolution induced by 
the Hamiltonian Hq on the system of landmarks x — (xi, • • • From a mechanical 

point of view, the landmarks represent the positions of n particles in and each is 
the momentum attached to the particle Xi. During the evolution, the particles interact 
and the equation = — (p, x) implies that the time derivative of the momentum of 

the i^^ particle is equal to the internal forces — ^^(p, x) acting on it. This is a simple 
extension of the usual Newton equations. 

Remark 1. Note that when there is no interaction between the particles i.e. Ky{z^z') = 
Id^d, then Ho{p,x) = ^Yl^^i and the system ^ reduces to mxi = pi and 
Pi = mXi = so that the particles evolve independently along straight lines at constant 
speed. Though this kernel is not admissible, it is not even continuous, this case can be view 
as a limit case when the characteristic scale X ^ in When the particles interact, 
the evolutions of the particles are no longer straight lines. The velocity of a given particle 
is the aggregate of the contributions coming from every particle as described in 

This Hamiltonian formulation can be efficiently retrieved by the application of the 
Pontryagin Maximum Principle (PMP)[lJ. The previous assumption in ([5| implies that 
the following control problem is controllable (since for every v = ('^i, . . . , '^n) ^ ^^^^ 
there exists u G V such that v = u ■ x = {u{xi))i<i<n)- For any chosen (initial and 
final) configurations of landmarks, there exists a path between the two landmarks sets. 
Following an optimal control point of view, we can introduce the Hamiltonian of the 
system 

(8) H{p,x,u) = {p,x) - ^{u,u)v = {p,u-x) - ^{u,u)v 
minimizing H w.r.t. to u we obtain: 

(9) u{.) = k{.,x)p, 

with the compact notation /c(., x)p = Yl7=i ^vi-^ Xi)pi. Then the minimized Hamiltonian 
can be written as: 

^ n ^ 

(10) Ho{p,x) = - ^{pj,Kv{xj,Xi)pi) = -{x,k{q,q)x) . 

Hence we get the Hamiltonian equations of the previous theorem by directly applying the 
PMP. 

It is rather standard to prove that this approach endows the landmark space (set of 
groups of distinct points) with a Riemannian metric which is induced from that of the 
group of diffeomorphisms. More generally once a positive definite kernel is given, it gives 
rise to a Riemannian metric on the space of landmarks. We then use the term kernel 
metric to designate such a Riemannian metric on the space of landmarks. 

The Hamiltonian framework (also viewed from an optimal control viewpoint) can be 
generalized to shape spaces where the landmark space L is replaced with with L^{R/Z, R^) 
for closed curves in the plane [13] or measures. It has been also generalized to discontinu- 
ous images [34J. However when deforming embedded objects in through the action of 
smooth vector fields, exact matching usually impossible. The associated control problem 
is not controlable any more. This is the reason why the Hamiltonian equations are estab- 
lished for the inexact matching problem which includes a penalty term in the minimization. 
The effect of this penalty term is often to smooth the structure of the momentum used in 
the Hamiltonian formulation. 



3. Spline interpolation on landmark space: a generative growth model 

3.1. The standard growth model. The usual growth model formulation as derived in 
[27] was initially defined in the case of images. Assume that one has a continuous time 
sequence of noisy image data if ^ t e [0, 1], and an image template /q. The basic growth 
estimation problem is to estimate a path (j)f, t ^ [0, 1] in the diffeomorphic group Gy (a 
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sequence of deformations) such that 1^ (j)t ■ Iq with (j) ■ I = I o (j) ^ . From a Bayesian 
perspective, this inverse problem is cast as the minimization of the cost 



(11) J'{v) = \j'^ \vt\ldt \IP - rt ■ h?dt 

where is constrained to be the flow of v starting at Id at time and v G ^^([0, 1], V^). 
Note that the underlying stochastic model for v is white noise in L2([0, 1], V) so that the 
associated flow is a Kunita flow [23j with almost nowhere different iable trajectories 



In the landmark setting addressed in this paper, assume that we have a sequence x 



D 

t 

of data (where = (xf^J is an n-tuple of d dimensional landmarks) and a template 
xq. a straightforward adaptation of the above growth model leads to the new cost to be 
minimized 

1 \ 



(12) ^'(^) = 2 J, X 

where (j) - xq = (0(xo,i)), and the squared distance \xP — (j)^ ■ xo\'^ comes from a stan- 
dard Gaussian white noise model the on individual landmarks. As above the underlying 
stochastic growth model is a Kunita flow for (j)t with irregular trajectories (see Fig. [t]). 
However, considering that x^ is continuously observed in time, the minimization prob- 



lem (12) can be cast into an optimal control problem associated to the reduced Hamilton- 



ian H given by [24J 

H{p,x,t) =Ho{p,x) - ^|xf -x\^ where Ho{p,x) = ^^pj K{xi,Xj)pj . 

The associated optimal trajectories are solutions of the Hamiltonian flow: 



(13) 



.Mo I 



The optimal time dependent vector field v can be reconstructed through the equality 

■^t(^) = ^K{z,xt,i)pt,i. 

i 

When A ^ 0, the optimal evolution converges to a geodesic evolution given by the 



usual Hamiltonian Hq. This is not surprising since in (13) the second term of the right 
hand side vanishes and leads to the minimization of the kinetic energy. 

Note that the filtered trajectory x, obtained from the observation x^, is in time 
for continuous time observations. It thus has the advantage of correcting the poor prior 
stochastic growth model given by Kunita flows. However, this is due to the fact that we are 
assuming continuous time observations. In contrast, for discrete observation times t/^, the 
optimal solution is piecewise geodesic, again no more C^, and offers limited interpolation 
properties. A better prior on x is needed. 

3.2. Perturbative growth model. Denoting Ut = \{xt — x^) the perturbation from 
the geodesic evolution, the estimated trajectory x extracted from the noisy observation 
xf*, t G [0, 1] can be reconstructed from the integration of 

r xt = ^{pt,xt) 

(14) I 

derived from the ODE ([T3| if one knows the initial momentum po and u. The fundamental 
idea is that this (ODE) can play the role of a generative engine for trajectories if we 
put the proper constraint on u seen as a control variable. Interestingly, from a mechanical 
perspective, Ut can be interpreted as external forces acting on the landmark configu- 
ration {ut^i is the force acting on particle i). In absence of external force, the landmark 
configuration follows a geodesic evolution. When the external forces do not vanishing, the 
trajectory followed by these landmarks deviates from a simple geodesic to various interpo- 
lation trajectories. The question of the dynamics of these external forces is of importance. 
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The first step is to note ttiat wiien Ut = X{x — xf)^ tiie cost of tiie filtered trajectory x 
can be written as 

(15) .r{v) = \j'^ \vt\ldt + ^ 1^*1'^*' 

SO tiie minimization of provides a control on both terms of the right hand side, both 
positive quantities. Moreover, for a given distribution P^d on x^, the associated distri- 



bution Px on the filtered trajectories x can be defined through the generative model (13) 
for the adequate distribution Pq on (po, u). Very little can be said on the distribution Pq, 
but as a first approximation we can proceed as follows If the overall model minimization 



of (15) for any observation x^ provides a reasonable filtered trajectory x and since for 
such X, the quantity ^ \ut\^dt < J^{v) should stay small (at least controlled), one can 
put a constraint on its expectation under Pq 

(16) Eo = Ep,{f \ut\^dt). 



Using no more information and maximizing entropy fiO\ |2Q] under the constraint (16), a 
reasonable first order model for the marginal distribution of u is given by a standard 
white noise i.e. Uf = adBt where it is reasonable to set clos^to the value A. 

The marginal distribution on po is more problematic, as well as the conditional distri- 
bution of u given po. The simplest solution is to assume independence of po and u: 

Po{dpo, du) = Po{dpo) (8) Po{du) 

so that that conditional generative model {po fixed) for the trajectory x can be defined as 
the SDE 

= ^{puxt)dt 



(17) 



dpt = -^{puxt)dt + adBt 



3.3. New growth estimation minimization problem. This new growth model gener- 
ates solutions and can be the foundation for a new formulation of the growth estimation 
problem in the realistic situation of sparse discrete observation times. Indeed, if we con- 
sider M observation times ti, • • • , ^m, this new prior can be used to derive a new growth 
estimation method : 

M 

|2 



(18) 



miJ{p^,u) = \E{po)^ \ [ \ut\'dt^^J2K 



k= 



subject to (x,p) solution of the ODE ( p3| ) 

where (po, u) G x i^^([0, 1],R^'^) with initial conditions (xo,Po) and E comes from the 
log-likelihood of the prior Po{dpo) on po. 

Remark 1. The choice of the regularization term E for po does not seem important in 
the finite dimensional case we are considering here. An improper flat prior can be used 
(as we do in the application below) and the term can be dropped. Alternatively consider 
Po fixed to 0. The weight on u is however of particular importance since it controls the 
deviation from a geodesic evolution corresponding to u = 0. 



f 

Jo 



One can consider more general penalties on u such as 

{Kx^ut,ut)dt 

/o 

for a state dependent metric {K^ is assumed here to be a positive symmetric matrix). 
However, there should be some rationale for the final choice. 

At this point we should mention the Riemannian cubic spline point of view developed 
in [291. In this framework, we start from a finite dimensional Riemannian manifold M 



The different approximations done in this derivation should encourage us to be a Httle cautious about 
any strong statement! 
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and the basic problem is to interpolate between two configurations {xq^xq) and (xi^xi) 
with a smooth curve 7 minimizing an extended bending energy : 

min^ 5(7) = \V^7?^dt 
subject to 

7(0) = xo, 7(0) = ^0, 7(1) = ^1 and 7(1) = xi 

where V is the Levi-Civita connexion and | I7 is the metric given by the manifold at 
the current location 7. It is quite clear that in the situation M = with the flat 
Euclidean metric, we get 5 (7) = \^\'^dt corresponding to the classical cubic splines 
[31] |2]. Moreover, one can check (see appendix) that in our Hamiltonian framework, we 
have 

(19) u=p^d,Ho = K-^\/:tx 

where = {Ky{xi^Xj))ij and is the metric tensor at the current landmark posi- 
tions. Thus, the bending energy is given in our situation as a function of u : 

B{u)= I {K^u,u)dt= I \u\l^^dt 
Jo Jo 

where \u\x,^ is the induced metric on the cotangent space. 

Remark 2. This penalization is different from the simpler LP' norm we derived previ- 
ously. On one hand, the Riemannian cubic splines are completely defined by a unique 
metric underlying the Riemannian structure of the manifold, keeping a pure geometric 
and intrinsic point of view. However, this intrinsic point of view is not as natural as it 
may look at first glance since it links tightly the internal forces (given by the term dxHo) 
and the external forces ( given by u) to the same metric structure. This is quite question- 
able since the internal (resp. the external) forces proceed from intrinsic (resp. extrinsic) 
phenomenons. 

We advocated in the previous subsection to link the choice of the metric on the control 
u to the noise model on the observation data. We derived the particular case of the 
white noise situation, which produces the standard Euclidean flat metric on the control. 
However, more general situations could be of some interest leading to non-standard metrics 
on u. 



3.4. A non standard metric. Assume that there exist a Hilbert space W and a smooth 
mapping : W^^ VK* such that the data term in ( 15 ) is replaced by 



where | \w* is the dual norm on the dual space W*. This situation is quite natural if we 
consider that the true observed quantity is not x^ but an element iIj{x^) G W*. This 
setting is of particular importance in the situation of shape modelling where the individual 
label of the different landmarks cannot be observed and the point- wise correspondence be- 
tween the template and the observation is problematic. The so-called measure framework 
developed in [12j makes intensive use of such a data term where ip{x) = - Xir=i is the 
empirical distribution of the landmarks and W is a proper Reproducible Kernel Hilbert 
Space. 

In this new case, the associated optimal trajectories are solutions of the Hamiltonian 
flow : 



(20) 



Pt = -^(Pu^t) + A(^'(xt))t(^(xf ) - ^{xt)) 



where ip'{xt)'' : W* K"'' is the Hilbertian adjoint of the differential ^p'ixt) : W'^ 
W* of at X. As above, denoting Ut = X{'tp'{xt))^ip{xP) — ip{xt)) and assuming that 
ip' (xt)^ ip' (xt) : R*"^ M"*^ is invertible (i.e. ip'{xt) is one to one ) along the trajectory Xt, 
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we get 

A 



1 



r) I |-r\™r / TK-^-o/wv — 



= ^j\w{xt)H'{xt)]-^UuUt)dt, 



where tt^^ is the orthogonal projection in W* on Im(?/^'(x)) which is the tangent space at 
of the sub- manifold of W* defined by the immersion ?/^. In particular, we can use 
as local metric = [ip' {x)"^ ijj' {x)]~'^ which is now position dependent and defines a new 
Riemannian metric on the Landmark space. 

In the mentioned situation of measure representation of landmarks we have the following 
proposition : 

Proposition 1. Assume that the kernel Kw associated with the RKHS W is . Then 
ijj : W^^ W* defined by tlj{x) = ^Yl7^i^Xi is differentiahle and [iIj' {x)"^ ijj' {x)] is the 
block diagonal matrix : 

1 d^Ky 
dxjdxj ' 



Proof. See appendix. □ 

4. Existence results and Euler-Lagrange equation for shape spline 

estimation 

In this section we provide several rigorous results concerning the shape spline frame- 
work. Indeed, although the overall picture is formally quite clear, it is necessary to provide 
rigorous statements about the main objects we have just introduced. In particular, con- 



ditions for the existence of solutions, global in time, of the perturbed evolution (14) need 



to be specified, as well as an existence theorem for the growth estimation minimization 



problem (18). On a more practical side, the minimization problem will be solved by gra- 
dient descent, for which we will compute the directional derivatives of functional J and 
the associated Euler-Lagrange equations. 

4.1. Existence of controlled evolution and weak dependency in the control 
variable. First define the function f{q^u) = {dpHo{x^p)^ —dxHo{x^p) -\- u)^ for q = 
(x^p) G X and u G W^^. We will consider the following hypothesis on V and its 
kernel Ky: 

HO : V is continuously embedded in Co(M^,]R^) and its kernel Ky is in each of 
its variable. 

Following is the first result on the existence in time of a solution to the perturbed evolution 



(14). 



Proposition 2 (Existence of solutions, global in time, of the controlled system). Assume 
(RO). Then for any u G L^([0, T], R^^) and for any initial condition qq = (xo,Po) ^ 
^nd X ^nd^ ^/^g^g g^^^^^ ^ ^ ^ C([0,T],R^'^ X R^^) such that 

^^^^ I xt = xo ^ J^dpHo{xs,Ps)ds 

Moreover^ there exists Ci , C2 , C3 > independent of u and qo such that for any t <T we 
have 

(1) HoixuPt) <Ht = Ci{Ho{xo,po) ^Tj^ \us\^ds), 

(2) \xt\ < \xo\^C2THy\ 

(3) \pt\ < (IpoI + /o^ \us\ds)eMC:^THy^). 

Proof. Under (HO), we have existence and uniqueness locally in time of the solution to the 
system ( [21] ) for any initial condition q^ = {xo^po). The only point to be checked is that the 
solution does not go to infinity in finite time. Let to < T be such that we have a solution 



qt = {xt.Pt) defined on [0,toLFrom ([21]), we get that Ht = Ho{xt,Pt) = Hq{xq,Pq) 

^3 = 



jQ{dpHo{xs,Ps),Us)ds. From (b|, since dp^Ho{xs,Ps) = v{xi) for v{.) = Ej^i ^v(., ^j)ft, 
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we get \dpHo{xsjPs)\ ^ CHo{xs^PsY^'^ for a universal constant C and we deduce that 
Ht < Ho{xo,po) + C(max5^[o,t] i^s)^/^ /q Hence 

(22) m^xHs<2Ho{xo,Po)^^C\f \us\dsf <2Hq{x^,p^) ^ AC^T j \usfds<HT 

for Ci = max(2, 4C^). The upper bound does not depend on to and in particular the 
Hamiltonian can not explode in finite time. It is sufficient now to prove that Xt and pt 
stay also bounded. Indeed, we have 

(23) \xt\ < \xo\ + I / dpHo{xs,Ps)ds\ < \xo\ f H^^ds < \xo\ + CTH^^^ . 

Moreover, \pt\ < \dxHo{xs,Ps)\oo\Ps\ds^ \us\ds-^\po\ and since again \dx^Ho{x,p)\ = 
\dv'^{xi)pi \ for 'u(.) = Y^^=i ^v{--, ^j)Pj-> we get from (|3| that there exists C > such that 
\pt\ ^ C' /q Hl^'^\ps\ds + Jq \us\ds + \pq\. Using Gronwall's Lemma, we get eventually 



\pt\ < (bol + |w.|ds)exp(C' f Hl'^ds) < (IpoI + K\ds)exp{C'TH^/^^ 
Jo Jo Jo 



(24) 

From ( [22] ), ( [23] ) and ( [24| ) we deduce that Xt and pt stay uniformly bounded on [0,to[. 
Since to < T is arbitrary, the system (21) admits a solution on [0,T] satisfying points 1), 
2) and 3) of Proposition [5] □ 

Proposition 3 (Dependence in u). Assume (HO). For any qo G M^^ x W^^ and any 

u e L2([0,T],R^'^), let G C([0,T],R^^ x M^^) (denote t/ie solution of ^ with initial 
condition qq. Then the mapping {qo-,u) q^'^^ is continuous for the weak topology on 
L2([0,T],R^^) and uniform convergence on C([0,T],R^'^ x R^^). 

Proof. Let be a weakly converging sequence in LP' and q^ q^^ be a converging 

sequence of initial conditions. There exists R > such that sup^>o \u^\2 ^ R and 
by Proposition [2j if q^ (resp. q"^) denotes the solution of (21) for (qo^u) = {qQ^u^) 
(resp. {qo,u) = {qo^.u'^)), then there exists M > such that Igfl < M for (t, n) G 
[0, T] X N U {oo}. Let Km > such that 

\f{q,u)-f{q\u^)\<KM\q-q'\^\u-u'\ 

for Ig''! < M and li, li' G R"^"^ (such Km exists since (HO) implies that dHo is C^). We 
have 

\q? = \f fiq:, <) - fiqT.Ods + {q^ - q^)\ 



< Km f - qT\ds + | A< - uT)ds\ + \q^ - q^\ 
Jo Jo 

so that using Gronwall's Lemma we get 

(25) kr - < (ko -q^\ + sup I I\k - uf )ds\)eMKMT) . 

r<T Jo 

Since u Ugds is a continuous linear mapping, we get from the weak convergence 

that I - uf)ds\ ^ as n ^ oo. Noticing that | - uf)ds\ < 2yjr' - rR for 

any r < r^ < Ascoli's Theorem turns the previous simple convergence into uniform 
convergence in r. □ 

Theorem 2 (Existence of a minimizer in the inexact case). Assume (RO) and that 
q ^ Xq e Mc^(R) is with {XqU^u) > c\u\'^ for some fixed c > and that C{q^u) = 
\{%qU^ u) . Assume that qo E{qo) is a non negative and lower semi- continuous function 
such that E +oo when \qo\ oo, that gk is a continuous non negative function for any 
1 < k < M and let xq G R^'^ x R^^ be a fixed initial condition. Then the function 

k=i 

defined for any {qo,u) G R^^ X L2([0,T],R^^) (where (7^°'^ is the solution of fm with 
initial condition qo ) reaches its minimum. 
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Proof. Let {Qq ^u'^)n>o be a minimizing sequence for J. Since the ^^'s are non negative 
and E > with E +00 when \qo\ +00, we get that {%q^u^ .,u^)ds and I^qI 
are upper bounded. Since we assume that (XqU^u) > c\u\'^, we deduce that is a 
bounded sequence in and we get by weak compactness of the strong bahs that (up 
to the extraction of a sub-sequence) there exists u"^ such that u'^ u'^ . Again, up to 
the extraction of a subsequence, we can assume that there exists q^^ such that Qq q^^. 
Since 

Jo Jo 
Jo 

we deduce from Proposition [3] that 

sup llXgr. -Xgoo|| ^ and / {{Xqr. - Xq^)u'l .u'Dds ^ . 

s<T ^ ^ Jo 

Moreover, by weak convergence {Xq'^u'f' ^ — u'f')ds ^ so that C{q'^^ u'^)ds < 

hminf C{qg^ u^)ds and by continuity of the QkS and of E^ J{Qq^^ u'^) < hm J(go 7 = 
inf J. □ 

Obviously, this proof gives also the existence of a minimizer in the exact case where the 
spline is constrained to go through a sequence {x^^)i<k<M of landmarks configurations if 
there exists at least one such controlled path with finite cost. 



4.2. Directional derivatives and Euler-Lagrange Equation. 
Proposition 4. Assume (EQ) and let u,6u e L^([0, T], R"'^) and qo.^qo ^ 



For any e > 0, we denote qt^e ^ C'([0,T], 
Ut,e = Ut -\- eSut and the initial condition qo^e = <7o + eSqo. Then we have 



the solution of (21) for the control 



(26) 



lim sup 

t<T 
e^O - 



^qt 







where dq is the absolutely continuous solution on [0, T] of the linearized system 

(27) Sqt = dqf{qt, Ut)Sqt + duf{qt, Ut)Sut 
with initial condition Sqo. 

Proof Assume < e < 1 so that supo<g<i < 00. As 

we did in the proof of Proposition [s] we deduce from Proposition |2] that qt^e = {xt,e^Pt,e) 
is uniformly bounded by some M > for (t, e) G [0,T]x]0, 1]. Again, if Km > is such 
that 

\f{q,u) - f{q',u')\ < KmIq - q'\ + \u - u'\ 

for l^l, \q'\<M and u, u' € K"-^ then 1^^,, -qt\<KM /q \qs,e - qs\ds + e{dqo + /q \Sus\ds) 
and by Gronwall's Lemma 

(28) \qt,, - qt\ < eexp{KMT){dqo + [ \Su\sds) . 

Jo 

Since from (HO), / is and dqf{q,u) and 9„/(g, u) are uniformly bounded for \q\ < M, 
there exists a unique solution 5q of (27) which is absolutely continuous. Moreover, 



I ft. 



qt 



(29) 
(30) 



Sqt\ 

f 

Jo 

f 

Jo 



f{qs,e,Us,e) - f{qs,Us) 



dqf{qs,Us)Sqs + duf{qs,Us)Su, 



ds 



dqf{qs,Us){ 



Qs,e 



ds 



f 

Jo 
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where 



OqJ[Qs,Us) h duJ[qs,Us)du, 



■ dqf{qs,Us)— 



However, from (28) and the fact that dqf is uniformly bounded for \q\ < M, we get that 
rjs^e is uniformly bounded on [0,T]x]0, 1]. Since rjs^e ^ for e ^ 0, we get by Lebesgue's 
Dominated Convergence Theorem that T]s,eds 0. Using (29) and Gronwall's Lemma, 
we get 

qt,e - qt 



Sqt 



< 



( / r]s,eds)exp{ / dqf{qs,Us)ds) 0. 

JQ Jo 



□ 



Theorem 3 (Directional derivative). Assume (RO), assume that q %q G M(^(R) is , 
that Qk is for any 1 < k < M and that E is . Let C{q^u) = ^{XqU^u) and T > tM- 
Then for any u,Su e ^^([0, T], R^^) and qo.Sqo e M^^ x R^'^, if 

nT M 

J(e) = £;((7o,e) + / C{qt,e^ut,e)dt^y29k{qtk,e) 



have 



(31) 



lim 



J(e) - J(0) 



M 



= (VE.Sqo)^ / {dqC{qs,Us)Sqs ^ duC{qs,Us)Sus) ds ^^{Vgk{qtk)i^qtk) 

(32) = {VE{qo)^Po,Sqo)^ [ {V uC{qs,Us) ^ duf^ {qs.Us)Ps,5us)ds 

Jo 

where C{q^u) = {%qU^u), Sqt is solution of (27) and Pt is of bounded variations with 
Pt = and 

M 

(33) dPt = -dqf{quUtfPtdt -Y,^9k{qt,) St, 

k=l 

where v <^ Sx denotes a vectorial Dirac measure at location x with value v. 

Proof For any variation 6u of the control we get by Proposition [4] that e qt^e is 
differentiable and deq\e=o = ^q where Sq = dqf{qs^Us)Sqs + duf{qs^Us)Sus. Moreover, 

Aie) = ^^M^ 



- {\/E{qo), Sqo) + / {dqC{qs,Us)Sqs + duC{qs, Us)Susds + ^(V^, 
V k=i 

- / dqC{qs,Us) 
Jo 



k{qtk)^^qtf,) 



M 



X^(V^/c(gtJ, 



qs,e - qs 



Qtk,e - qtk 



Sqs ds 



k=l 



f 

Jo 



where 



M 



C = E f ^^^^^^^^^^^ - {VgM, ^^^^^)V^^^^^^^^^^-(Vi5(*),^9o) 



and 



' dqC[qs,Us)^-^ ^ duC[qs,Us)ou, 



From the fact that the g^s and E are and (28), we get that Ce ^ with e ^ 0. Since 
qs^e is uniformly bounded for (s, e) G [0,T]x]0, 1] and %q is C^, one easily gets from (28) 
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that there exist a, 6 > such that \r]s^e\ < <^ + ^(I'^ls + l^'^sP)- Since r]s^e ^ as e ^ 0, 
we get by Lebesgue's Dominated Convergence Theorem that rjs^eds — ^ 0. Using (26), 
we get eventually 74(e) so that 

dJ , 



(34) 



SJ=^{0) = {\/E{qo),Sqo) 



nT M 

+ / {dqC{qs,Us)5qs^duC{qs,Us)5us)ds^^^{Vgk 

"^0 k=l 



and (31 ) is proved. Introducing now Mt^s the semi-group solution of dgMt^s = dqf{qs, Us)Mt^s 
with Mt^t ^^2nd we get 



and from (34) 



(35) 



Thus we have 



5J 



Sqs = / Mt^s dufiqt, Ut)Sutdt + Mo^s^qo • 

^0 



6J={\/E{qo),Sqo)^ / {\/uC{qt,Ut),Sut)dt 



k=l 



{VuC{quUt)^duf{quUtf / MT^V qC{qs,Us) ds 



where 



+ duf{qt, ut)^ Mt^tJ'^gk{qtk)U<tk , ^ut)dt 
k=i 

+ (V^(^o) + / Ml^VqC{qs, Us) ds^y2 ^lt,^9k{qt,)U<t,.hQ) 

k=i 

= [ {^uC{qu ut) + duf{qu UtfPt, Sut)dt + (V^(^o) + Po, Sqo) . 
Jo 

Pt = / Mj^^^\/qC{qs,Us)ds^^M^^,^\/gk{qt,)lt<t,- 

k=i 

One easily sees that Pf is absolutely continuous between the observation times with jumps 
at the observation times. Moreover, differentiating t -> Mt^t'Mf^t = Id, we get dtMt^t' = 
-Mt^t'dqf{qt,ut) so that we get ([33|. □ 

From Thm|3j we get immediately the Euler-Lagrange equations for shape splines evo- 
lution : 

qt = f{qt,ut) 



(36) 



dPt = -dqf{qu UtfPtdt - Y.k=i 5. 



tk 



where Pt = (Pf Pf) with boundary conditions 



(37) 



Pt = 0, VE{qo) + Pq = . 



Remark 1. Note that in our experiments, we will consider that xq is fixed and let po 
be free (sometimes called natural spline in the classical cubic splines framework). This 
corresponds to a flat prior on po and gives the new boundary conditions 

(38) PT = o,pP = o. 

Moreover, the cost functions g^ will usually not depend on p so that dpgk{q) = and Pf 
is absolutely continuous in time. In particular, we get in this case that for natural splines, 
Ut is continuous and vanishes at the boundary of the interval [0,T]. 
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5. Numerical experiments 

In this section we provide preliminary experiments illustrating the behavior of shape 
splines in simple 2D synthetic experiments. We start with a small family of observation 



1.5 




-1.5 -1 -0.5 0.5 1 1.5 



Figure 1. Left hand side :2D plots of the sequence of sampled curves 
Ck for k = 0, • • • ,5. The initial curve is a circle and the final shape is 
a horizontally pinched ellipse (plain curve). The intermediate curves are 
generated by linear interpolation between the initial and final shapes. 
Note that no noise is added here and that the intermediate shapes are 
rescaled by a factor depending on k. Right hand side : estimated 
shape spline displayed in space-time representation. The color indicates 
the value of the norm of u on the surface through time (blue for low 
values, red for high values) 

times < ti < • • • < ^ ^ and for each time, a landmark configuration G M^"^ 
defined as a noisy version of regularly sampled curves c^. To stay close to some realistic 
framework where the observation points are sparse, the number M of time points is kept 
small to M = 5 (with = 3.8). In these experiments, we focus on the simple Euclidean 
metric — p in the data term, which is in agreement with our Gaussian noise 
assumption. 

In our first experiment, we consider the evolution of an initial 2D circular shape which 
somewhat linearly evolves into a pinched ellipsoidale shape. The shapes are regularity 
sampled in a consistent way so that the problem of point correspondence between two time 
points does not need to be consider (see Fig[T]). To emphasize the space-time regularity 
provided by the spline shapes interpolation, we display the result as a surface in the 3D 
space-time. The first interesting point to note is that the shape spline provides an actual 
smooth interpolation of the evolution in time between the observation epochs and extends 
to shape spaces the specific behavior of classical cubic splines. In particular, the shape 
spline actually joins the observed shapes with very good accuracy despite the fact that 
we are using here inexact matching. A second important point is that the shape spline 
comes with the estimation in time of the control variable u which can be interpreted as 
an external force bending the underlying geodesic. What we see in this example and in 
the other similar situations, displayed in Fig [2j is that the point-wise values of u give 
interesting information on the evolution process. More specifically, in every example, the 
initial and end shape have specific features (numbers of lobes, orientation, etc) and the 
most active zones are easily interpretable and correspond to transition regions in the shape 
evolution. 

5.1. Robustness to noise. The robustness to noise is a rather important subject from 
a practical point of view. Indeed, the noise basically degrades the spatial resolution 
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Figure 2. The figure displays two estimation experiments of shape 
splines for two different problems. From left to right, the first column 
displays the sequence of observed shapes, the second and the third col- 
umn two different viewpoints. The color of any point on the space-time 
surface is related to the norm of its control variable (with increasing val- 
ues from deep blue to red) 

of the measurements so that the evolution through time of a particular point of the 
evolving curve may be a sharply broken line. The standard spline approach can be quite 
efficient in filtering this noise if the time sampling frequency is high enough. This is hardly 
the case in many important situations. However, neighboring points behave coherently 
through time and offer an interesting source of spatial redundancy. Much of the large 




Figure 3. Robustness to noise. An i.i.d. Gaussian noise with standard 
deviation <j = 0.1 is added to each measurement point. Shape splines 
are computed with two different scale parameters, A = 0.001 on the left, 
A = 0.6 on the right. 

deformation shape space theory involves the integration of spatial redundancy in the 
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comparison between shapes. The shape sphne setting, considering shapes as a whole 
and not as a bag of independent points, keeps this important aspect but adds a time 
component and considers the problem in the full space-time setting. This robustness to 
noise is illustrated in Figure |3] where an i.i.d. Gaussian noise with standard deviation 
cr = 0.1 is added and a series of shape splines are computed under increasing values of the 
spatial regularity scale parameter A as introduced in Q . For low values of this parameters 
(with respect to the overall scale of the shapes) the reconstructed evolution is clearly far 
from any reasonable solution since the spatial redundancy is hardly taken into account. 
Increasing the value of A to values in accordance to the scale of the object produces a 
much better reconstruction of the actual shapes at any observation time but also keeps 
existing time regularity. 

5.2. Extrapolation. Another distinguished feature of the usual spline setting which is 
extended in the shape spline setting is the fact that the extrapolation of the data outside 
the interval of observation is quite straightforward. Indeed, outside the limits of the 



2 




Figure 4. Extrapolations. The extrapolation of the evolution at both 
ends of the observation interval (A = 0.6, a = 0). 



observation interval, the value of the control parameter u is set to zero and the evolution 
is naturally extended with a geodesic evolution. Moreover, one can check (see Remark 
[l]) that u vanishes at the last observation time so that the previous extension is for 
the control variable u and for the shape variable x. Note that in the standard growth 



model described in (12), the evolution is extrapolated by fixing the shape variable to its 
value at the last observation and this extrapolation is only . In Figure [ij we display 
a simple example of extrapolation where the underlying evolution within the observation 
interval is the same than in Figure [T] The computed extrapolation appears visually quite 
natural at both ends. 

5.3. Comparison with piecewise geodesic evolution. We end this section with a 



comparison with the piecewise geodesic interpolation scheme [27J, derived from (12) for a 
finite set of observation points as the minimizing solution of 

M 

(39) j''{v) = ^ I \vt\ut+iy\xr,-xt. 



= 1 f \vt\Ut+if2K 

-^^ i=i 



where as previously stated, Xt = (j^ti^o) and (j)^ is the flow of G L^([0, 1], y). We display 
an experiment on a synthetic evolution where t x^ is given by a simple analytic formula 
where a circle evolves smoothly into an ellipse by increasing its eccentricity through time 
combined with a rotation of the principal axis (see three orthogonal views of the synthetic 
object in the middle column of Fig. [5|. We display the estimated evolution in the piecewise 
geodesic setting given by (39) and with a shape spline in FigjsJ As expected, the piecewise 



geodesic estimation provides good results but with a loss of regularity at the observation 
points as in the simpler situation of piecewise linear interpolation in signal processing. 
The shape spline seems to perform better at the observation points but also to provide 
a better estimation between observation points. In Fig. |6j we provide a more quantified 
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Figure 5. Comparison with piecewise geodesic evolution. On this ex- 
ample the middle column corresponds to 3 different views of a synthetic 
shape evolution (here the time axis is vertical in the first two rows and 
horizontal in the last row). The first column correspond to the piecewise 
geodesic evolution one can get from (12) and the last one to the shape 
spline estimation. 



10 




piecewise geodesic 
shape spline 



10 



10 




10 



Figure 6. Comparison of the error E (see (40)) between piecewise 
geodesic and shape spline estimation. The horizontal axis is the number 
M of observation points and the vertical axis the distance in arbitrary 
units through time between the estimated shape and the actual shape 
provided by the synthetic evolution. Right panel, log- log plot. 



comparison of the approximation quality between the estimation process by computing 
the error 

(40) \J -^tfd'^ 

as a function of the number M of observations. It is quite clear that the convergence is 
faster with shape spline interpolation than with the piecewise geodesic interpolation (we 
do not go beyond 11 observation points since the error is small enough to be approaching 
other numerical errors in the optimization scheme). The loglog-plot in Fig. [g] seems 
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to indicate a polynomial convergence in C jM^ very similar to the classical situation in 
interpolation theory that (a = 2 for linear spline and a = 4 for cubic spline 



6. A STOCHASTIC SHAPE SPLINE MODEL 



In this section, we study further the first candidate for the second order model (17) of 
stochastic evolutions of shapes. We prove that the solutions are well defined for all times 
and present some simulations that highlight some features of this model. 

6.1. Non blow-up result. This section is devoted to study the well-posedness of the 



SDE (17) introduced as our generative growth model in subsection 3.2 The random force 
is chosen to be an increment of the Brownian motion though we could also have introduced 
a Levy process in the evolution of the momentum to account for sudden activations of 
cells. It would turn our growth model into a more realistic one. However, this Brownian 
perturbation is the first step toward such a model and we will now discuss its feasibility 
from the mathematical point of view. We will prove that the solutions of the SDE do not 
blow up in finite time a.s. 

The stochastic differential system is: 

(41a) dpi = -dxHo{pt, Xt) dt + edBt , 

(41b) dxt = dpHo{pt, Xt) dt . 

Here, £ is a constant parameter and Bt is a Brownian motion on R^^. We will work 
with the Gaussian kernel but this can be directly extended to other kernels. When e 
is constant, there is no difference studying these stochastic differential equations with 
the Ito integral or Stratonovich one. However, for a general variance term, we will use 
the Ito stochastic integral. From the theorem of existence and uniqueness of solution of 
stochastic differential equation under the linear growth conditions, we can work on the 
solutions of such equations for a large range of kernels. Yet in our case the Hamiltonian is 
quadratic, and the classical results for existence and uniqueness of stochastic differential 
equations only prove that the solution is locally defined. In the deterministic case, this 
quadratic property could imply existence of a blow-up. To prove that the solutions do 
not blow up in finite time in the deterministic case (e = 0), we can use the fact that the 
Hamiltonian of the system is constant in time. By adapting that proof (also closely related 
to the proposition [2| and controlling the Hamiltonian, we will prove that the solutions 
are defined for all time. 

A first remark we will use is the following, for any a G and z G M^, 

(42) {a,Kv{z,z)a)j^d < C^Ir^ • 

Thus, we introduce the stopping times defined as follows: let M > be a constant and 

(43) TM = {t > I max(|x,|, \pt\) > M} , 

let also Too = limM^oo t be the explosion time. Differentiating Ho{ptATM ^ ^tATM) with 
respect to t, we get on {t < tm)'- 

dHo{t) = dxHo{pt,Xt)dxt dpHo{pt,Xt)dpt -^^tT{Kv{xi{t),Xi{t))) — dt. 

i=l 

In the deterministic case the Hamiltonian is constant, whereas here the stochastic pertur- 
bation gives 

dxHo{pt, xt)dxt + dpHo{pt, xt)dpt = edpHo{pt, Xt)dBt . 

Thus 



/ dHo{t) = / e{dpHo{pt,Xt),dBt)^ / V tr(i^y (x,(t), x,(t)))-dt , 
Jo Jo Jo -^^ ^ 



and 

p2 

(44) E\R^{j)TArM.^TArM)\ < Ho{0) ^ E{C^ -dfiT A Tm) < i/o(0) + CVdnT. 



^Note that we have implemented here a least square approximation algorithm (see and \39\ and 
not an exact interpolation algorithm but with a value 7 weighting the data attachment term high enough 
to make the data error negligible. 
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Now, we aim at controlling Xiatm using the control on dxt given by \dpHo{pt^ Xt)\oo < 
C^/Ho{pt,Xt): 

(45) \Xr^^t\ < \xo\ + / CHo{Ps,Xs)^^^ds < \xo\ + / CHoip 

s/\tm 1 '^S/\Tm 

Jo Jo 

Jo 

However, < At P a.s. and by monotone convergence theorem (recall that Hq is non 
negative), 

/ rtATM \ 
E{At) = lim (|xo| +^ / CHo{PsArM.XsArM)^^^ds . 

Also, 

(46) e(^J^ ^ HoiPsATM^^sATMy^^dS^ -^{lo ^ob.ATM^^.ATM)'/'^^) 

Fub. J E (^HoiPsATM^^sATMy^^) ds < J E{Ho{Psatm^^sAtm))^^'^ ds 

CS+ (|44| 

We deduce 



< ViU (i^o(O) + C'^e^dns) ds 



1/2 



E{At) < \xq\ + CVi (^j^ (i^o(O) + C'^e^dns) ds^ < oo and < oo P a.s. 



and as a consequence 



lim sup \xtATM I ^ +^ P 

M^oo 



We also control the evolution equation of the momentum as follows, 

rtATM rtATM 

(47) \ptATM \ < / \dMPs.Xs)\ds + Ipo + / sdBs 

Jo Jo 

Now we use the assumption (|3| to control dxHo{p^x): 

\dMP,x)\ < \p\\dv{x)\ < C\p\h'o^\ 



We rewrite inequality (|47|) and we use Gronwall's Lemma to get: 

ptATM rtATM 

\PtATM\< / C\ps\Ho{ps,Xsy^'^ ds^\po^ / sdBs\, 
Jo Jo 

IPM.M I < (\Po\ + sup I sdB, \) efr- CHoir>.,^.V^^^s 

\ u<t Jo J 

bM.„|<fbo|+ sup I /^dB,|)er^-<^^o(P--)^^^'^^ 

\ U<tAToa Jo / 

The first term on the right-hand side \po\ + sup^<^/\^^ | edBg \ is bounded by \po\ 
sup^<^ I edBs I < ooPa.s. and with inequality (46) we have that 

Since on (too < t) one has 

lim maxd^MTMl. \PtATM\) = Ji"^ = oo, 
we deduce P(roo ^ t) = and Too = +oo almost surely. 
We have proved for the case e{p^ x) = eld^ 
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Theorem 4. Under assumption the solutions of the stochastic differential equation 
defined by 

dpt = -dxHo{pt,Xt)dt^e{pt,Xt)dBt 

dxt = dpHo{pt,Xt)dt. 

are non exploding when £ : R^^ x R^^ ^ L(R^^) is a Lipschitz and bounded map. 

Proof To extend the proof to the case when e is a Lipschitz and bounded map of p and 
X, we can prove that the preceding inequahties are stih vahd. 

First, with the Lipschitz property of e the solutions are stih defined locahy. The Ito 
formula is now written as, on (t < tm) 

dHo{t) = dxHo{pt, xt)dxt + dpHo{pt, Xt)dpt + ^tr(£^(p^, Xt)Kx,e{pt, Xt))dt . 

where is block matrix defined by = {Ky{xi,Xj))i<ij<n' 
We still have the inequality ( [44| ) with 

tT{e^{pt,Xt)K:,^e{pt,Xt)) < {Cnd\e\oof 

if \s{p,x)w\'^ < \s\1o\w\'^ where |e|oo denotes the supremum norm. Indeed, if {ei)i^[i^n(i] 
the canonical basis of R^^, denoting s = s{x^p)^ we have 

nd nd 

i=i i=i 
where is the largest eigenvalue of K^. We have < tT{Kx) = tr(i^y (x^, x^)) 

and using (42) we get tr(i^y (x^, x^)) < <^AJ^^^^^ < dC'^ so that AJ^^ < C'^nd. Hence, 

nd 

tT{s'K,e) < C^ndY,\^\lc < {Cnd\e\^f . 

i=l 

Thus we get. 



/ 

Jo 



dHoit) < / {dpHo{puXt),e{puXt)dBt) + / ^ ^^dt, 

Jo Jo ^ 



E[Ho{T A tm)] < ffo(O) + ^( (^^^Hoo)" TAtm)< Ho{0) + {Cnd\e\^f T , 

and all the remaining inequalities follow easily thanks to the control on Hq and the bound 
on e. □ 

Once this stochastic model is well-posed on landmark space, the question of its exten- 
sion to shape spaces naturally arises. It can be proved that this stochastic model does 
have an extension to the infinite dimensional case: in the case of L2(R/Z,R2), the natural 
extension of the Brownian motion on the landmark space is a cylindrical Brownian motion 
on L^(R/Z,R^). It may be somewhat surprising to deal with such irregular noise on the 
momentum variable, we stress the fact that this noise is read by the kernel which strongly 
regularizes the noise. Though we will not develop it further, it proves that this model has 
a consistent extension to continuous-shape spaces and could be used to deal with random 
evolutions of continuous shapes. 

6.2. Simulations. With these simulations we illustrate the interesting features we ob- 
served above. First, this model gives realistic perturbations of geodesies contrary to a 
first order model such as a Kunita flow. A simulation of a Kunita flow is illustrated in 
figure Fig. [7] where the evolution of 40 points on the unit circle is represented under a 
Gaussian kernel of width 0.9. The time evolution has, as expected, the roughness of a 
Brownian motion and the space variation is smoother due to the kernel. In comparison, 
our stochastic model gives smoother evolutions in time as in Fig. [To] and Fig. [ll] 

Second, our stochastic model is a perturbation of a geodesic evolution and this nice 
property is illustrated in figures Fig. [9] - \TT\ 

Figure Fig. [8] shows the geodesic evolution of 40 equidistributed points on the unit 
circle for a Gaussian kernel of width 1.0, the target configuration for the landmarks is 
obtained through a simple affine transformation that gives the final ellipse. On these 
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Figure 7. A simulation of Kunita flow with 40 points on the unit circle on the left 
of the flgure. The z axis (blue arrow) represents the time. 

simulations the color change only represents time. Figures Fig. [9| - pT] represent stochastic 
perturbations of the previous geodesic; we progressively increase the standard deviation e 
of the noise from ^Jne = 0.9 to 1.7 and finally ^/ne = 3.5 (the noise is rescaled w.r.t. the 
number of landmarks to converge to a well defined SPDE at the limit (see [34J)). Each of 
these three figures represents one Monte- Carlo simulation of the stochastic model with a 
simple Euler scheme. 



Figure 8. Geodesic evolution 
- White unit circle as initial 
shape. 



Figured, white noise pertur- 
bation of the geodesic (same ini- 
tial momentum po), y/ne = 0.9 



Figure 10. increasing the FIGURE 11. Increasing the 

variance of the noise, ^/ne = 1.7 variance of the noise, y/ne = 3.5 



The simulations in figures Fig. 12 show the position at time 1 of the 40 points for 5 
Monte-Carlo simulations. On the two figures we plotted the initial momentum po (at- 
tached to the 40 points) associated with the geodesic from the initial circle to the target 
ellipse. As this model was designed to produce random shape evolutions, it can also be 
used as a generative engine to produce random shapes. Increasing the noise also increases 
the expectation of the energy of the system since the Ito formula applied on the Hamil- 



tonian in subsection 6.1 shows a linear growth in time of Hq proportional to e. This can 
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be guessed when comparing the two displayed cases in Fig. [12] since in the second one 
the noise is 4 times bigger. For any statistical estimation of the model parameters, this 
property should be somehow taken into account. 




Figure 12. Effects of the standard deviation e of the noise : 5 simulations of 
random deformations of the unit circle. The initial momentum is fixed, the kernel 
width A = 1 (see Left-hand side : y^e = 0.25 ; right-hand side : y^e = 1. 

An important feature of this model is that the noise is "read" by the kernel. We show 
in Fig. [13] simulations of the model for a null initial momentum on the same initial shape 
and we decrease the width of the Gaussian kernel from 3 to 0.3. The standard deviation 
of the noise is constant set to ^/ne = 1.0. 




Figure 13. Effects of the width of the Gaussian kernel. Left-hand side :A = 0.3 ; 
right-hand side : A = 3.0. 

These last simulations show the importance of the choice of the kernel and as a by- 
product the choice of the operator e in front of the noise will be also important. Now we 
can formulate a stochastic model for evolutions of shape that would be closer to realistic 
evolutions: 

^^^^ \dpt = -dqHo{pt,qt) -^ut^ edBt 

where Uf is of bounded variations. At this point, we underline that the model param- 
eterization is completely open and it should be tightly related to consistent statistical 
estimations. 
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7. Shape splines on homogeneous space 

In this section we provide a more formal and geometrical picture of what could be an 
extension of the shape spline to the previously mentioned important cases. For this, we 
need to introduce some of the standard vocabulary of geometrical mechanics as developed 
in [25]. We will try as much as possible to avoid the conceptual burden of the intrinsic dif- 
ferential and symplectic geometry through the extensive use of local coordinates. Readers 
looking for a more intrinsic formulation could refer to [25j. 

7.1. Geometrical setting. The proposed framework for shape spline is given by three 
ingredients: a group G of transformations, a Riemannian manifold Q and a left action of 
G on Q denoted (^, q) Lq{g) = g ■ q. We will assume also that for any q ^ g ^ g ■ q 
is a surjective submersion (i.e. the differential of Lq has full rank everywhere). 

7.1.1. Local coordinates. Basically G will be a (finite dimensional) Lie group with Lie 
algebra (5 on which we consider a right invariant metric given by a dot product on (5. The 
infinite dimensional setting where G is a group of diffeomorphisms is more involved and 
requires more analytical work as in [33]. This is clearly the target setting we have in mind 
but we want in this paper to stay away from any complicated analytical developments. 
To keep the focus on the global picture, we will assume implicitly that we work in a finite 
dimensional setting for which the existence of all the introduced objects is straightforward. 

Denoting q = (^^, • • • ^q^) local coordinates on Q and since {dq^^ • • • ^dq^) is a basis 
of T*(5, we can write any a G T^Q as a = ^Pidq^ so that (g^, • • • • • • ,Pn) are 

local coordinates on the cotangent bundle T*Q. Given q = (g^, • • • ,^^), we will denote 
P = {Pir " iPn) a generic element of T*Q and m = (q^p) a generic element of T*Q as we 
did previously in the fiat case of landmarks. 

7.1.2. Infinitesimal actions and cotangent lift. The first thing we need is to extend the 
action of G on Q to an action of G on the cotangent bundle T*(5. Note that the differen- 
tiation in g of g ^ q ■ q ai g = Idc yields an infinitesimal action {£,^q) ^ ^ • g for ^ G 0. 
Differentiation in q yields the action (^, Sq) g ■ Sq G Tg.qQ for Sq G TqQ and by duality 
the action (g^p) ^ g ■ p ^ "^g gQ P ^ ^qQ^ uniquely defined through the equality 

(49) {g-p\g-Sq) = {p\Sq). 
We denote 

(50) {g,m) ^ m- g = {g ■ q,g ■ p) 

for m = {q,p) G T*Q the induced action on T*Q. 

In summary, the initial action g ^ g • q on Q is naturally lifted to an action g ^ g • m 
on the cotangent space T*(5 (usually called cotangent lift [25j). Differentiating the action 
(g^m) ^ g ■ m on the cotangent space T^Q ai g = Mg, we get the infinitesimal action on 
T*(5, which is defined in local coordinates by ^ • m = (^ • ^ • p) where 

(51) {^■p\Sq) + {p\^-Sq)=0 



as obtained by differentiation of the conservation equation (49) 



7.2. Euler-Poincare equation. The initial matching problem between shapes in Q is 
defined as the solution of the optimal control problem with fixed boundary 

™n^t h Io{Lit^it)dt 
subject to 

Qt = ' Qt^ QO — ^init, Ql = ^targ 

where L : 6* ^ 6 is the isometry between (3 and its dual induced by the metric on 0. 

The Hamiltonian associated to the classical matching problem is given in local coordi- 
nates by H{q,p^^) = {p\^ • q) — ^(L^|^) with reduced form 

H{q,p) = ^{KJ{q,p)\J{q,p)) 
where J{q^p) G 0* is uniquely defined by 
(52) {J{q,Pm) = {p\^-q) 
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for any ^ G (usually called the momentum map). To obtain the Hamiltonian evolution, 
we need to compute the variation of 5 J as a function of the variation Sq and 5p in q and 
p. Introducing in local coordinates the so-called symplectic matrix 



Idn 
-Idn 



one checks easily that {SJ\^) = • Sq) + {Sp\^ • q) so that using (51) we get 

(53) {5J\0 = -(C • P\5q) + {Sp\^ ■ q) = -(J(C • m)\5m) 

and dH = —J{KJ • m). Since the associated Hamiltonian evolution equation are given by 



m 



IdH we get 



(54) 

or equivalently 
(55) 



m 



Kj • m with j = J{q,p) 



q = i - q and p = (, - p with ^ = KJ{q,p) . 



This extremely simple expression of rfi in term of the momentum map and the infinitesimal 
action on the cotangent space reveals part of the nice geometrical structure underlying 
the evolution. In this setting, it is interesting to consider the time evolution of the pair 
(g, j) instead of the pair {q^p) since j follows an autonomous equation. Indeed, from ([51]) 



and (|55|), we get that for any C G we have (f | C) = (6 • P I C • ^) + {P I C.(6 ' q)) 



{p I Cj^t • ^) - 6 • (C • q)) = -{P I ad^,(C) • q) = -{j \ ad^,(C)) where ad^(C) 
adjoint representation of the Lie algebra (3 so that we get 



K,C] is the 



dt 



ad 



KjJ 



0. 



This equation, called the Euler-Poincare equation plays a central role in geometric me- 
chanics and more recently in the large deformation methods in shape analysis and com- 
putational anatomy [TBI 123 • 

Now consider the perturbed dynamic p = Kj{m).p + u with u G T^Q or equivalently 

-+ad^,, = /. 

where h = J{q,u) and the associated optimal control problem for the state variables 
(QiJ) G Q X 6* and the cost ll^^lg^^. The norm we consider here is the dual norm 

induced on T^Q by the metric on TqQ. Let %q : T^Q ^ TqQ he the isometry such that 
{u I %qu) = \u\q^^. The Hamiltonian associated to our new control problem for the costate 
variable {pq, () eT^Q x (3 is given by 

(56) :K{qJ,PqX^u) = {pq I Kj ■ q) + (-ad^^^j ^ h \ () - ^{XqU \ u) . 



To compute its reduced form let us note that from (51), we get -^^ih \ () = ( ■ q. Hence 



du 



'K = implies XqU = C ' Q giving the reduced Hamiltonian 



(57) 
(58) 



nqd,p,, = (P, I Kj ■ q) + {-sd*Kjj I C) + 2 IC • q\l 

= {J{q,p,) I Kj) + (i I ad^i^i)) + \\C-q\l, 



where | \q denote the norm on TqQ given by the Riemmanian metric on Q. The associated 
Hamiltonian evolution is derived quite easily: 



(59) 



dq 
dt 
dpq 
dt 



dt 
dC 
dt 



Kj 
= ^-q 

- ^-Pq = 

-ad^(j) : 
Fadc(C) 



icC2\C-q\l) 



KJiq,p,)=0 



ifadj(i) 

Here again, the derivation should be considered at a formal level or in a smooth finite 
dimensional setting since existence results are beyond the scope of this paper. 
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8. Conclusion 

In this paper, we present new tools for shape evolution analysis or growth analysis 
in computational anatomy through the introduction of second order evolutions. Shape 
splines seem to overcome some of the limitations of the previous first order schemes and 
with their stochastic counterpart can provide the backbone of new statistical time regres- 
sion tools. From various perspectives, shape splines offer new interesting mathematical 
and practical challenges: extensions to the infinite dimensional case of continuous shapes 
and to images, development of efficient and scalable numerical schemes to solve the spline 
estimation problem, integration of time realignment as developed in [9j, derivation of more 
complex stochastic engines beyond the white-noise situation presented here and develop- 
ment of consistent statistical schemes in the spirit of [3j. The solutions to some of these 
problems appear well within reach. 

In this paper we preferred to stick to the finite dimensional setting in order to focus 
on the general picture and avoid more difficult analysis, the more technically involved 
infinite dimensional situation has been partially explored in the stochastic case in [34J. 
Acknowledgments. The authors would like to thank Darryl D. Holm and Colin J. 



Cotter for a fruitful discussion about introducing noise on the momentum variable which 
was the starting point of the stochastic model presented here. 

9. Appendix 

9.1. Link with cubic splines. We recall that on a Riemannian manifold (M, ^m), a 
cubic spline between (xq^vq) and (xi^vi) is a curve c : [0, 1] ^ M that minimizes 



1 

2 J gM{'^cC,Vcc)dt . 



(60) J(c) 2 

The Euler-Lagrange equation for this functional is the following (see [29l fTT|) 

(61) V?c + i?(V^c,c)c = 0. 

Let c : / ^ M a smooth path, then the covariant derivative is given in coordinates by 

VcC= + ^C/er*^^^(c)Q]ie[l,r], 
k,l 

and we observe that the second member of the right hand side only depends on (c, c). We 
compare this expression with the Hamiltonian equations (here M = £): 

(62) c = ^k{c,c)p = [dik{c,c)]{c)p^k{c,c)p, 

(63) c = [dik{c^ c)](/c(c, c)~^p)p — k{c^ c)dxH + A:(c, c)u , 
where for any x^y G R^^, k{x,y) denotes the block matrix defined by 

k{x,y) = iKv{xi,yj))i<ij<n- 
Now the geodesic equations are given by 

Ci^'^CkTi i{c)ci = Oforz G [l,r], 

k,l 

where ^ = gjzi^ didj^dk) are the Christoffel symbols in the chosen coordinates. 

However, if ?x = we obtain the geodesic equations in the Hamiltonian form. Then we 
can identify 

-\^(^k'^k,M(^i\ie[i,r] = [dik{c,c)\{k{c,c)~^p)p-k{c,c)dxH , 

k,l 

which gives 

VcC = /c(c, c)u . 

Now we have 

gL {"^cC, Vcc) = gL {k{c, c)u, k{c, c)u) = {u, k{c, c)u) , 
which proves the desired result. 
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9.2. Proof of Proposition [l| Let us denote 9i2^vf(^i, ^2) = ^^^^(^15^2) for any 
^1,^2 G and Az^h,e = — ^2) for any G and e > 0. Let us recall that for 



any ^1,^2 G 



(64) {5^^ , ) VF* = (^1 ^ ^2 ) . 

First, we start with the proof that Az^h,e converges in when e ^ 0. Indeed we get 
from (l64l that 



(65) (A^^,^^,ei, A^2,/i2,e2)H^* = di2Kw{zi^ seihi,Z2-i-te2h2) - hi^h2dsdt 

Jo Jo 

(66) = di2Kw{zuZ2) ■hi^h2^ o{{\eihi\ + |e2/i2|)|/ii||/^2|) , 

so that \Az^h,e — Az^h,e'\w* ~ ^(kl + W\)' Since 1^* is complete, Az^h,e converges in 
to a limit point denoted (5^^^ such that z^^hi^ ^2^2^^* = lim(A^i,/,,,e, A^2,/i2,e)Ty2 = 
9i2^vf(^i, ^2) - hi <S) h2. In particular ^ ^ is a continuous linear mapping from W^^ 
to ly*. 

Now, we get from (65) 

\d.,h-{S. + Si^^)\l,=o{\hf) 
SO that i/j is differentiable at any location x G M.^^ with differential ?/;'(x) given by 

1 



Thus, using again ( 65 ), we get {iIj' {x)Ax,iIj' {x)Ax)w* = 7^ Xli^j=i ^12^(^*7 ^jO'^^^^^^i 
which proves Proposition [T] 
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